Technotools/Technotools (Chestnut CD-ROM)(1993).ISO/lang_c/cug188/trans.doc

Text File | 1985-08-21 | 14KB | 270 lines

.m⌠ 0
.mΓ 11
.p∩ 13
/*
HEADER:        ;
TITLE:        C elementary transcendentals;
VERSION:    1.0;

DESCRIPTION║ááá    Sourcσááá codσááá fo≥áá al∞áá standarΣááá ├                            ì
               transcendentals

        Employ≤á ldexp(⌐á anΣá frexp(⌐á functions╗á iµ                ì
               suitablσá version≤á oµ thesσ arσá no⌠á provideΣ                ì
               b∙á ß giveε compiler¼á thσ version≤ provideΣ iε                ì
               sourcσá codσá wil∞ requirσá adaptatioεá t∩á thσ                ì
               doublσ floa⌠ format≤ oµ thσ compiler.

KEYWORDS:    Transcendentals, library;
SYSTEM:        CP/M-80, V3.1;
FILENAME:    TRANS.C;
WARNINGS║áá    frexp(⌐áá anΣá ldexp(⌐á arσáá implementatioε                ì
               dependent«á  Thσá compile≥á employeΣá doe≤á no⌠ ì
               suppor⌠ááá minu≤áá (-⌐áá unar∙áá operator≤áá iε                ì
               initialize≥á lists¼á whicΦ arσ requireΣ b∙á thσ                ì
               code.
AUTHORS:    Tim Prince;
COMPILERS:    MIX C, v2.0.1;
*/
              Transcendenta∞ Functioε Algorithms

     Al∞á programminτ language≤ whicΦ arσ descendeΣ froφá FORì
TRA╬á o≥á Algo∞ includσ pre-defineΣ transcendenta∞á functions«     ì
StandarΣá textbook≤ oε programminτ assumσ tha⌠ thσ usua∞á matΦ ì
function≤á arσ alread∙ availablσ anΣ no⌠ wortΦ thσá pape≥á reì
quireΣá t∩ describσ ho≈ the∙ migh⌠ bσ coded«á  Thi≤ i≤á unforì
tunate¼á a≤ fe≈ compute≥ system≤ providσ reliablσ transcendenì
tals¼á anΣá ß basiπ understandinτ oµ thei≥ content≤ i≤ helpfu∞ ì
iε thei≥ use.

     Thσá followinτá example≤ usσ polynomia∞á representations¼ ì
becausσá the∙ leaΣ t∩ compac⌠ codσ anΣ givσá averagσá accurac∙ ì
withiε onσ bi⌠ oµ thσ bes⌠ available«  Faste≥ algorithm≤ exis⌠ ì
bu⌠ requirσ morσ codσ o≥ datß storage¼á iµ writteε iε thσ samσ ì
language«á  Commercia∞á version≤ oµ thσ function≤ arσá usuall∙ ì
writteεá iεá assembler¼á bu⌠ thσ idea≤ arσ bes⌠á expresseΣá iε ì
highe≥á leve∞ languages«á  Obscurσ trick≤ whicΦá aεá anonymou≤ ì
code≥ though⌠ woulΣ improvσ speeΣ arσ useles≤ wheε thσ prograφ ì
fail≤á anΣ thσ sourcσ codσ i≤ no⌠ available«á  Thesσ function≤ ì
havσ beeε testeΣ iε nearl∙ thσ samσ forφ oε severa∞á computer≤ ì
iε botΦ FORTRA╬ anΣ C.

     Thσá author'≤ hopσ i≤ tha⌠ thesσ example≤ ma∙ se⌠ ß miniì
muφ standarΣ fo≥ accurac∙ anΣ speeΣ oµ ├ librar∙á transcendenì
tals«á  Failinτá this¼á thσá individua∞ use≥ shoulΣá havσá thσ ì
opportunit∙á t∩ checδ whethe≥ an∙ deficiencie≤ oµ hi≤ codσ ma∙ ì
bσ duσ t∩ thσ librar∙ provideΣ witΦ thσ ├ compiler«á C'≤á sooε ì
t∩á bσ standardizeΣ suppor⌠ oµ bi⌠ field≤ anΣ lo≈ leve∞ operaì
tion≤ oε double≤ help≤ iε thσ expressioε oµ thesσá algorithms«  ìèSincσá C¼á eveε iε thσ futurσ ANS╔ standard¼á wil∞ no⌠ suppor⌠ ì
purσ singlσ precisioε oε computer≤ whicΦ werσ no⌠ designeΣ fo≥ ì
mixeΣ precisioε arithmetic¼ anΣ thσ standarΣ librar∙ function≤ ì
arσ double¼ wσ sho≈ onl∙ thσ doublσ precision¼ witΦ thσ underì
standinτ tha⌠ change≤ iε thσ polynomial≤ woulΣ adap⌠ t∩á othe≥ ì
precisions«á  Thσ filσ TRANSLIB.FO╥ show≤ FORTRA╬ version≤á iε ì
singlσ precision¼ whicΦ shoulΣ servσ a≤ ß demonstratioε oµ thσ ì
change≤ requireΣ t∩ var∙ precisioε anΣ language.

     Thσá polynomia∞á coefficient≤ werσ determineΣ b∙á runninτ ì
onσá oµá Stevσ Ruzinsky'≤ program≤ (1)«á  Iεá eacΦá case¼á thσ ì
polynomia∞á i≤ deriveΣ b∙ truncatinτ ß standarΣ Taylo≥á serie≤ ì
o≥á rationa∞á polynomia∞ representatioε oµá thσá functioεá anΣ ì
usinτ Steve'≤ prograφ t∩ obtaiε morσ accurac∙ witΦ fewe≥ term≤ ì
iε thσ requireΣ range«á  Iε somσ cases¼ 2╡ digi⌠ arithmetiπ i≤ ì
requireΣá t∩ fi⌠ ß polynomia∞ accuratσ t∩ a⌠ leas⌠ 1╢á digits«  ì
Table≤ oµ coefficient≤ havσ beeε generateΣ fo≥ polynomial≤á oµ ì
an∙ accurac∙ requireΣ u≡ t∩ 2░ o≥ morσ significan⌠ digits.

     Trigonometriπá function≤ caε bσ calculateΣ witΦá adequatσ ì
efficienc∙á usinτ portablσ code¼á eveε iε archaiπ dialect≤á oµ ì
FORTRAN¼á Pascal¼á anΣ C«á  Wheε usinτ polynomials¼ i⌠ appear≤ ì
bes⌠á t∩ usσ thσ sin(⌐ functioε a≤ thσ basiπ serie≤ anΣ calcuì
latσ thσ other≤ froφ it«á  ┴ polynomia∞ fo≥ tan(⌐ require≤á a≤ ì
man∙ term≤ a≤ sin(⌐ anΣ cos(⌐ together«  Sincσ cos(⌐ anΣ tan(⌐ ì
ma∙ bσ calculateΣ froφ sin()¼á onl∙ onσ shorte≥ tablσ oµ coefì
ficient≤á i≤á needed«á  Iε effect¼á tan(⌐ i≤ calculateΣá b∙á ß ì
rationa∞ polynomia∞ approximation«á  Expressinτ thesσ trigonoì
metriπ relationship≤ b∙ #definσ teache≤ thσ preprocesso≥ algeì
brßá whicΦá aεá optimizinτ compile≥ coulΣá usσá t∩á advantage«  ì
Checδá tha⌠ you≥ compile≥ doe≤ no⌠ executσ function≤á multiplσ ì
time≤á a≤ ß resul⌠ oµ macr∩ expansion«á  Tan(⌐ coulΣ bσ calcuì
lateΣ mucΦ faste≥ witΦ ß biτ tablσ a≤ i≤ buil⌠ iε t∩ thσ 8087¼ ì
bu⌠á mos⌠ program≤ wσ havσ seeε spenΣ morσ timσ oεá sin(⌐á anΣ ì
cos().

     Thσá firs⌠á ste≡ iε calculatinτ sin(⌐ i≤á t∩á reducσá thσ ì
rangσ t∩ |xⁿ ╝ pi/▓ b∙ subtractinτ thσ neares⌠ multiplσ oµ pi«  ì
Thσá ODD(⌐ functioε showε assume≤ tha⌠ (int⌐ wil∞ extrac⌠á thσ ì
lo≈á orde≥á bit≤á iε casσ thσ argumen⌠á exceed≤á maxint«á  Thσ ì
multiplicatioεá anΣ subtractioε shoulΣ bσ donσ iε lonτá doubleô ì
iµ i⌠ i≤ available¼ bu⌠ iε an∙ casσ thσ rangσ reductioε shoulΣ ì
no⌠ reducσ accurac∙ wheε thσ origina∞ argumen⌠ wa≤ alread∙á iε ì
thσá desireΣá range«á  Thσá methoΣ mos⌠ ofteε useΣá save≤á onσ ì
dividσ (whicΦ ma∙ bσ changeΣ t∩ ß multiply)¼ a⌠ thσ cos⌠ oµ aε ì
unnecessar∙á roundoff«á  Worse¼á thσ resul⌠ ma∙á underflo≈á t∩ ì
zero«á  Underflow≤á iε thσ codσ showε occu≥ onl∙ wherσ ßá zer∩ ì
doe≤ no⌠ affec⌠ thσ result«  WitΦ thσ precaution≤ taken¼ extrß ì
precisioεáá i≤á no⌠á needeΣá anywherσá othe≥á thaεá thσá rangσ ì
reduction.

     Mos⌠á compute≥á system≤ includσ ßá polynomia∞á evaluatioε ì
function¼á sometime≤ iε microcode«  T∩ savσ spacσ wσ prefe≥ t∩ ì
usσ sucΦ ß functioε eveε thougΦ i⌠ ma∙ bσ slowe≥ thaεá writinτ ì
ou⌠á thσá polynomia∞ iε Horne≥ form«á  SucΦ function≤á operatσ ìèlikσ poly(⌐ showε here¼á bu⌠ (presumably⌐ faster«á  Thi≤ moduì
larit∙á permit≤á specia∞á hardwarσ feature≤ t∩á bσá useΣá morσ ì
effectively¼á althougΦ thσ timσ requireΣ t∩ cal∞ anothe≥ funcì
tioε ma∙ bσ excessive«  Loo≡ unrollinτ (showε here⌐ o≥ odΣ anΣ ì
eveεá summinτ (seσ exp2()⌐ arσ ofteεá effective«á  Thσá serie≤ ì
useΣá fo≥ elementar∙ function≤ havσ term≤ increasinτ iε magniì
tudσ fas⌠ enougΦ s∩ tha⌠ therσ i≤ n∩ neeΣ fo≥ extrßá precisioε ì
intermediatσ computatioε.

     Thσá arctaε functioε use≤ severa∞ substitution≤ t∩ reducσ ì
t∩ ß rangσ oµ 0 < x < 1«á  AlthougΦ thσ Taylo≥ serie≤ doe≤ no⌠ ì
convergσ ove≥ thi≤ range¼á thσ minima° polynomia∞ does¼ bu⌠ s∩ ì
man∙ term≤ arσ needeΣ tha⌠ onσ furthe≥ ste≡ oµ rangσ reductioε ì
i≤ preferred«á  Codσ lengtΦ ma∙ bσ tradeΣ fo≥ speeΣ b∙ usinτ ß ì
tablσ oµ tangent≤ witΦ ß shorte≥ series«á  ┴ rationa∞á polynoì
mial(3⌐á i≤ useΣ b∙ man∙ run-timσ libraries¼á a⌠ somσ cos⌠á oµ ì
accuracy«á  ┴á serie≤ witΦ onσ les≤ terφ woulΣ stil∞ givσ ove≥ ì
1╡á digit≤á accuracy¼á whicΦ i≤ thσ maximuφ possiblσá oεá somσ ì
systems.

     AlthougΦá loτ basσ ▓ i≤ no⌠ usuall∙ includeΣ iε thσá lis⌠ ì
oµá librar∙ functions¼á i⌠ i≤ almos⌠ alway≤ useΣ a≤ thσá basiπ ì
logarithφá fo≥ binar∙ floatinτ point¼á becausσ i⌠ i≤ thσá mos⌠ ì
efficien⌠á basσ fo≥ exponentiatioε witΦ ß floa⌠á power«á  Man∙ ì
implementation≤ fai∞ t∩ obtaiε thσ logarithφ oµ number≤á closσ ì
t∩ 1«á  Thσ followinτ routinσ i≤ thσ mos⌠ widel∙ testeΣ oµ thσ ì
example≤á giveε here¼á anΣ ha≤ beeε founΣ accuratσ oε al∞á maì
chine≤ whicΦ subtrac⌠ ▒ accuratel∙ (somσ don't!)«á  Iε theory¼ ì
precisioεá i≤á los⌠ wheε thσ exponen⌠ i≤ addeΣ iε a⌠ thσá end¼ ì
bu⌠á thi≤á i≤ no⌠ seriou≤ iε thσ mos⌠ commoε case≤ oµá number≤ ì
nea≥ 1«  Makinτ thσ functioε lonτ doublσ woulΣ correc⌠ this.

     Log▓á caεá bσ writteε iε portablσ codσá iµá thσá standarΣ ì
function≤á frexp(⌐ anΣ ldexp(⌐ fo≥ extractinτ anΣ changinτ thσ ì
exponen⌠ fielΣ oµ ß floatinτ poin⌠ numbe≥ arσá available«á  Oε ì
man∙á systems¼á adequatσ speeΣ wil∞ no⌠ bσ obtaineΣ unles≤ thσ ì
codσá i≤á supplieΣ iε line«á  Thi≤ shoulΣá bσá possiblσá usinτ ì
union≤ oµ structures¼á bu⌠ direc⌠ usσ oµ an∙ applicablσ floatì
inτ poin⌠ hardwarσ instruction≤ i≤ preferable«á  Thσ algorithφ ì
require≤ separatσ calculatioε oµ thσ loτ oµ ß numbe≥ nea≥ ▒ b∙ ì
ßá series¼á whicΦ i≤ addeΣ t∩ thσ scalσ (neares⌠ intege≥ powe≥ ì
oµ 2⌐ oµ thσ argument.

     Usinτá bi⌠á field≤á oµ structure≤ a≤á aεá alternativσá t∩ ì
frexp(⌐ anΣ ldexp(⌐ caε bσ tricky¼á bu⌠ wil∞ allo≈ ß shor⌠ in⌠ ì
comparisoε t∩ bσ substituteΣ fo≥ thσ floa⌠ comparison« Thσ DE├ ì
systeφá i≤ unusua∞ iε tha⌠ floatinτ poin⌠ number≤á arσá storeΣ ì
witΦá bytσ pair≤ reverseΣ whilσ integer≤ arσ storeΣ witΦá ful∞ ì
reversσ bytσ order«á  Thi≤ ma∙ no⌠ ruiε thσ bi⌠ fielΣ code¼ iµ ì
n∩á fielΣ crosse≤ ß 16-bi⌠ boundary«á  Thσ ├ functioεá frexp(⌐ ì
conflict≤á witΦá thσá suggesteΣ IEE┼ standarΣá whicΦá give≤á ß ì
resul⌠ jus⌠ greate≥ thaε ▒ rathe≥ thaε jus⌠ les≤ thaε 1.

     ┴á commoε erro≥ iε commercia∞ implementation≤ oµ log(⌐ i≤ ì
t∩ bypas≤ thσ conversioε oµ ß comparisoε int∩ aε incremen⌠á oµ ìèthσá scalσá b∙ subtractinτ anΣ addinτ sqrt(.5⌐ insteaΣá oµá 1«  ì
Thi≤á ha≤á neve≥ beeε seeε t∩ worδ well¼á a≤ sqrt(.5⌐á i≤á no⌠ ì
representeΣ exactly¼ anΣ al∞ significancσ oµ aε argumen⌠ whicΦ ì
i≤á nea≥á ▒á ma∙ bσ lost«á  I⌠ isn'⌠ eveε a≤ fas⌠á unles≤á thσ ì
machinσá i≤ onσ whicΦ i≤ bette≥ suiteΣ t∩á extendeΣá precisioε ì
floatinτá poin⌠á thaε intege≥ anΣ logica∞á operations«á  Sincσ ì
mos⌠ system≤ havσ immediatσ constan⌠ representation≤ oµ 1« anΣ ì
2.¼ thσ differencσ iε codσ lengtΦ i≤ small.

     Iεá ßá ful∞ IEE┼ standard(3⌐á system¼á log(⌐á woulΣá neeΣ ì
specia∞á codσ t∩ takσ carσ oµ log(0⌐ anΣ infinitie≤ anΣ assurσ ì
tha⌠ accurac∙ i≤ no⌠ los⌠ oε subnormals«á  Thσ codσ showε wil∞ ì
givσá ß system-dependen⌠ resul⌠ fo≥ log(0⌐ whicΦ i≤ sucΦá tha⌠ ì
exp(log(0)⌐ stil∞ return≤ 0.

     Thσá exponentia∞á i≤ thσ mos⌠ difficul⌠ oµá thσá standarΣ ì
function≤ t∩ calculatσ witΦ ful∞ accuracy«á  I⌠ coulΣ bσ writì
teεá iεá portablσá codσ iµ therσ werσ ß wa∙ t∩á ascertaiεá thσ ì
exponen⌠ rangσ a⌠ compilσ time╗á IEE┼ P75┤ dictate≤ ß slightl∙ ì
differen⌠á rangσ froφ olde≥ system≤ anΣ require≤ specia∞á casσ ì
handlinτ oµ ou⌠ oµ rangσ arguments«á  Iµ ldexp(⌐ take≤ carσ oµ ì
over- anΣ underflow¼ thesσ check≤ ma∙ bσ omitted«  Man∙ variaì
tion≤ ma∙ bσ found¼á bu⌠ the∙ havσ iε commoε ß rangσ reductioε ì
baseΣ oε splittinτ thσ argumen⌠ int∩ aε intege≥ anΣ ß fractioì
na∞ part¼ raisinτ ▓ t∩ thσ powe≥ oµ eacΦ part¼ anΣ (iε effect⌐ ì
multiplyinτ t∩ ge⌠ thσ result.

     Unlikσ thσ othe≥ functions¼á whicΦ havσ n∩ eveε term≤á iε ì
thei≥ Taylo≥ series¼á exponentiatioε ha≤ ß symmetr∙ whicΦá caε ì
bσá exploiteΣá onl∙ iε ß rationa∞ approximation«á  Sincσá
2^(-x)=1/2^x¼á thσ rationa∞ polynomia∞ take≤ thσ forφ
P(x)/P(-x)¼á wherσ thσ numerato≥ anΣ denominato≥ arσ identica∞ ì
excep⌠ fo≥ thσ sigε oµ thσ odΣ powe≥ terms«á  Onl∙ halµ oµ thσ ì
term≤á neeΣ t∩ bσ evaluated.

     Thσá coefficient≤ ma∙ bσ calculateΣ b∙ penci∞á anΣá pape≥ ì
algebra¼á convertinτá ß continueΣ fraction(4)ô int∩ ßá rationa∞ ì
approximation(5)«  Sincσ thσ equation≤ fo≥ thσ coefficient≤ t∩ ì
fi⌠ thσ approximatioε t∩ ß se⌠ oµ point≤ takσ thσ samσ forφ a≤ ì
fo≥á ß simplσ polynomial¼á thσ minima° fittinτ prograφ ma∙á bσ ì
useΣ t∩ reducσ thσ numbe≥ oµ term≤ required«  Thσ precisioε oµ ì
thσá fina∞ resul⌠ canno⌠ bσ a≤ grea⌠ a≤ thσ precisioεá oµá thσ ì
numerato≥á anΣ denominato≥ polynomials¼á sincσ change≤ oµá thσ ì
las⌠á bi⌠á iεá numerato≥ anΣ denominato≥ changσá thσá quotien⌠ ì
irregularl∙ bu⌠ b∙ aε averagσ oµ tw∩ bits«á  Usσ oµ table≤á t∩ ì
reducσá thσá numbe≥á oµá term≤ requireΣ woulΣá speeΣá u≡á thi≤ ì
function.

     Sinh(⌐á shoulΣá havσ it≤ owε polynomia∞ fo≥ |xⁿá ╝á ▒á t∩ ì
maintaiε accurac∙ anΣ savσ thσ timσ whicΦ woulΣ bσ requireΣ t∩ ì
calculatσá eveε powe≥ term≤ iε exp()«á  Tanh(⌐ shoulΣá includσ ì
codσá t∩á givσ +-▒ directl∙ a≤ ß resul⌠ fo≥á largσá arguments«  ì
Thσá solutioεá usinτ ß specia∞ exp-▒ functioε i≤ workablσá bu⌠ ì
slowe≥á anΣ incompatiblσ witΦ thσ rationa∞á approximatioεá fo≥ ì
exp().è
.cp 5
     Sqrt(⌐á i≤á bes⌠ codeΣ a≤ aε arithmetiπ operatioε oεá thσ ì
samσ leve∞ a≤ divide¼á s∩ tha⌠ i⌠ caε executσ faste≥ thaεá anΣ ì
a≤ accuratel∙ a≤ divide«á  Man∙ system≤ fai∞ t∩ d∩ thi≤ anΣ s∩ ì
ß typica∞ procedurσ i≤ shown«  ┴ minima° polynomia∞ i≤ useΣ t∩ ì
ge⌠á aεá approximatioε accuratσ t∩ ▓ o≥ │ digits¼á anΣá Newtoε ì
iteration≤á arσ performeΣ unti∞ maximuφ accurac∙á i≤á reached«  ì
I⌠ i≤ faste≥ t∩ perforφ thσ maximuφ numbe≥ oµ iteration≤ whicΦ ì
ma∙á bσá requireΣ rathe≥ thaε t∩ g∩ unti∞ thσá las⌠á iteratioε ì
make≤ n∩ significan⌠ change«  Iµ thσ arithmetiπ operation≤ arσ ì
properl∙á rounded¼á thσ resul⌠ wil∞ bσ withiε onσ bi⌠ oµá ful∞ ì
machinσ accuracy.

     Thσá codσ requireΣ t∩ implemen⌠ thσ algorithm≤ appear≤ t∩ ì
bσá consisten⌠á witΦ publisheΣ standard≤ fo≥ thσá ├á language«  ì
Somσ compile≥ designer≤ conside≥ tha⌠ thσ Referencσá Manual(6⌐ ì
definitioεá i≤ satisfieΣ b∙ acceptinτ bu⌠ ignorinτ minu≤ sign≤ ì
iεá initializers«á  Thσ codσ wa≤ testeΣ b∙ patchinτá thσá conì
stant≤á iεá thσá compile≥ outpu⌠ code«á  Thi≤ probleφá ma∙á bσ ì
circumventeΣ b∙ expandinτ thσ polynomial≤ iε line¼á whicΦ i≤ ß ì
favorablσá tradσ oµ lengtΦ fo≥ speeΣ oε somσá machines«á  Thi≤ ì
sor⌠á oµá probleφá doe≤ no⌠ occu≥ eveε iεá TRANS   DOCσσ!╪⌠ code.

     Test≤ sho≈ tha⌠ FORTRA╬ anΣ ├ version≤ oµ thesσ function≤ ì
ma∙á bσ improveΣ iε speeΣ o≥ codσ lengtΦ b∙ assemble≥á coding«  ì
Excep⌠ possibl∙ fo≥ exp()¼á thσ accurac∙ oµ thesσ function≤ i≤ ì
a≤á gooΣá a≤ thσ bes⌠ anΣ frequentl∙ bette≥á thaεá thσá store-ì
bough⌠á versions«á
.paèReferences

1.Ruzinsky¼á Steveεá A«á  "┴ Simplσ Minima°á Algorithm,óá Dr« ì
Dobb'≤ Journa∞ #9│ pp« 84-91¼ Jul∙ 1984

2.Moshier¼á Stepheε L«á  "Compute≥ Approximations,ó  Bytσ pp« ì
161-178, Apr. 1986

3.Cody¼á W«á J.¼ Karpinski¼ R«"┴ ProposeΣ Radix- anΣ ì
Word-length-independen⌠á StandarΣ fo≥ Floating-poin⌠á Arithmeì
tic,ó IEE┼ Micr∩ pp« 86-100¼ Augus⌠ 1984

4.Abramowitz¼á M.¼á Stegun¼á I«á  "Handbooδá oµá Mathematica∞ ì
Functions,ó Dove≥ 1968

5.Press¼á W.H« e⌠ al«  "Numerica∞ Recipes,ó Cambridgσ Univerì
sit∙ Pres≤ 1986

6.Ritchie¼ Denni≤ M«  "Thσ ├ Programminτ Languagσ - Referencσ ì
Manual,ó  iε UNI╪ Programmer'≤ Manual¼á v.2¼ Hol⌠ Rinehar⌠ amΣ ì
Winston¼ 1983